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SUMMARY 


The coordinate transformation procedure and solution technique described 
herein have been shown to be well suited for the calculation of the complete 
flow field surrounding various two-dimensional and axisymmetric bodies. This 
conclusion is supported by the good comparisons obtained between predicted 
values and experimental data for pressure coefficient distribution and heat- 
transfer distribution on a cylinder and for density distributions and shock 
standoff distance on spheres. The solutions for the flow over the Viking aero- 
shell and a Jupiter probe were obtained by simply adjusting five transformation 
constants without having to perform any adjustment on boundary conditions. The 
versatility of the technique in being able to model many different axisymmetric 
blunt body shapes and its ability to calculate both the forebody and wake flow 
over these shapes represent a significant advance in aerothermodynamic 
technology. 

An upper limit on Reynolds number on the order of 1000 is imposed because 
of an inability to resolve a thin boundary layer with a reasonable number of 
computational mesh points. A lower limit on Reynolds number on the order of 
100 is imposed so that continuum theory is applicable. An upper limit on Mach 
number approximately equal to 4.0 is imposed because of a tendency to calculate 
negative static enthalpies in the process of capturing the bow shock for large 
Mach numbers. These limitations are a consequence of the numerical method used 
and not of the coordinate system itself. A more refined numerical method could 
ease some of these restrictions. 


INTRODUCTION 

The calculation of the complete flow field surrounding axisymmetric or 
two-dimensional blunt bodies in a supersonic stream has been a goal of fluid 
dynamicists which has only recently been made practically attainable through 
advances in high-speed computers. Inviscid shock-layer solutions (refs. 1 
to 4), viscous shock-layer solutions (ref. 5), and solutions of the full 
Navier-Stokes equations (ref. 6) have been obtained for the shock layer ahead 
of various blunt configurations. Many of those solutions have been specialized 
for various flow regimes which range from high Reynolds number flow (ref. 7) to 
low Reynolds number flow of a viscous rarefied gas (ref. 8). 

Mathematical descriptions of the flow in the wake of a blunt body are more 
difficult to obtain. Boundary-layer approximations which are successfully used 
in the shock layer are not applicable to a large portion of the wake where vis- 
cous forces are dominant (ref. 9); however, early predictions of the wake flow 
used these approximations (refs. 10 to 13)* Weiss (ref. 14) improved on this 
approximation and reported temperature field calculations which compared more 
closely with experimental data than previous work. Allen (ref. 15) solved the 
full Navier-Stokes equations in the steady state for the laminar wake behind a 


two-dimensional rectangular block. That solution specified a boundary-layer 
profile under a uniform stream as an inflow condition. 

This ’’patchwork" approach of solving the complete flow field, that is, the 
matching of separate, specialized solutions for the different areas of the flow 
field, was noted by Scala and Gordon (ref. 16). As an alternative approach, 
they presented a solution to the complete time-dependent Navier-Stokes equa- 
tions for the transient supersonic flow around a right circular cylinder. 
Kitchens extended that approach to obtain the steady-state solution to the same 
problem (ref. 17). 

It should be noted here that a patchwork approach for obtaining the com- 
plete flow field is a useful technique. In fact, these specialized techniques 
almost always lead to a more efficient, faster running numerical solution 
because fewer computational points are needed to define the flow field and many 
terms can be omitted from the Navier-Stokes equations by using an order of mag- 
nitude type approach (as is done in deriving the boundary-layer equations) in 
certain areas of the flow field. However, the major disadvantage of the patch- 
work approach is that it is often difficult to determine where and how the var- 
ious specialized solutions are to be matched. 

Existing solutions to the Navier-Stokes equations for the entire flow 
field have been obtained for simple body shapes. (See ref. 18.) An accurate 
description of the flow field over complex body shapes is facilitated by the 
choice of an appropriate coordinate system. Gnoffo (ref. 19) outlined many of 
the advantages and disadvantages that are associated with the use of various 
coordinate systems which can be used as the basis for the numerical differenc- 
ing of the Navier-Stokes equations over complex axisymmetric or two-dimensional 
body shapes. Some important considerations are the ability of the coordinate 
system to concentrate mesh points near the body for resolving the boundary 
layer and near regions of sharp curvature to treat rapid expansions. Orthogo- 
nality of coordinate lines to the body simplifies the boundary-condition 
specification. 

A generalized orthogonal natural coordinate system was presented in ref- 
erence 19 which gives close analytic approximations to various complex body 
shapes. The body shape approximations avoid the problems of discontinuous 
slope or curvature which occur in many aerodynamic configurations and yet are 
able to closely approximate these shapes. The purpose of this study is to show 
that this coordinate system can in fact be used to obtain the complete flow 
field over complex body shapes. 

A computational technique has been developed which uses this coordinate 
system for obtaining the solution of the Navier-Stokes equations for the entire 
flow field. A computer program has been written to describe flow over two- 
dimensional body shapes or axisymmetric body shapes. Comparisons with experi- 
mental data have been made to verify the technique and to ascertain where prob- 
lems due to the nature of the flow or due to numerical instabilities can occur. 
Flow fields around the Viking aeroshell and a candidate configuration for a 
Jupiter probe are calculated and presented. The analysis has been restricted 
to supersonic flow of a perfect gas with free-stream Mach number less than 4.2 
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and Reynolds numbers of the order of 100 to 1000, for reasons that will be dis- 
cussed subsequently. 
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SYMBOLS 


transformation constants, nondimensionalized by Rfj* 
speed of sound, nondimensionalized by V * 

OO 

drag coefficient defined in appendix F 


pressure coefficient, = 


p * - p * 

^b oo 

1 O 

- p *V * 2 

2 OO OO 


specific heat at constant pressure, J/kg-K 

error term defined in equation (22); also exponential 

metric coefficient, nondimensionalized by %* 

heat-transf er coefficient, W/m^-K 

total enthalpy, nondimensionalized by 

static enthalpy, nondimensionalized by 

thermal conductivity, W/m-K 

body length, nondimensionalized by Rjg* 

intermediate body length, nondimensionalized by %* 

free-stream Mach number, = V */a* 

7 OO 

integer in transformation equation 
Knudsen number, = X^/Rjg* 

Prandtl number, = y*c p */k* 

Reynolds number based on nose radius of curvature, 5 *Rjj*/y* 

total number of mesh points in 0-direction 
total number of mesh points in r|-direction 
integer in transformation equation 


3 



iiiiiiiiiHiiiiimim i 


n* distance normal to body, m 

p pressure, nondimensionalized by P 00 *V co *^ 

q c * convective heat transfer, W/m 2 

R radius of curvature, nondimensionalized by Rjj* 

R* gas constant, J/kg-K 

r ,0 ,<}5 coordinates in transformed space 

S* maximum cross-sectional area, m 2 

s arc length, nondimensionalized by Rfj* 

T temperature, nondimensionalized by T * 

T^* free-stream temperature, K 

t time, nondimensionalized by %*/V * 

U component of velocity in x -direction, nondimensionalized by V * 

u component of velocity in 0-direction, nondimensionalized by V^* 

Y component of velocity in y-direction, nondimensionalized by V^* 

VJf free-stream velocity, m/s 

v component of velocity in r-direction, nondimensionalized by V^* 

x,y,z Cartesian coordinates, nondimensionalized by R^* 

Y metric coefficient, nondimensionalized by Rjj* 

Y may maximum body radius, nondimensionalized by Rfj* 

a^- thermal accommodation coefficient 

3 coordinate stretching parameter in equation (17) 

Y ratio of specific heats 

A s * shock-layer thickness, m 

e smoothing parameter in equations (23) 

ri stretched coordinate in computational plane defined in equation (17) 
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m ai 


■ 1 min 11 


X* mean free path, m 

U viscosity, nondimensionalized by 

V m * free-stream viscosity 

p density, nondimensionalized by P^* 

P^,* free-stream density, kg/m3 

t shear stress 

4^ angle defining velocity vectors 

Superscripts : 
n time step index 

* signifies dimensional quantity (MKS system) 

Subscripts : 

aw adiabatic wall 

B base 

b body 

calc calculated 

exp experimental 

i index on mesh point in 6-direction 

j index on mesh point in n-direction 

N nose 

n index on transformation constants 

o refers to slip boundary condition 

r refers to conditions along constant 6 and 4> 

s shock 

stag stagnation 

6 refers to conditions along constant r and 4> 
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<j> refers to conditions along constant r and 9 

OO free stream 

Bars over symbols generally denote average values; however, in appendix B 
the bars denote variables at a predicted time step. (See step 1 of appendix B.) 
Arrows over a symbol denote vectors. 


COORDINATE SYSTEM 


Coordinate Transformation 


A generalized curvilinear orthogonal coordinate system which can be used 
for approximating various axisymmetric and two-dimensional body shapes is pre- 
sented in reference 19. The body shapes include spheres, ellipses, spherically 
capped cones, flat-faced cylinders with rounded corners, circular disks, and 
planetary probe vehicles. The transformation from the (0,r,4>) domain to the 
(x,y,z) domain for an axisymmetric coordinate system is written as 

N 

x(0,r,4>) = (-B sinh r + C cosh r) cos 0 - ^ A n e nr cos n0 

n=2 


y(6 ,r,4>) 


N 

(B cosh r - C sinh r) sin 0 + A n e nr sin n0 

n=2 


cos 4> 




( 1 ) 


z(0 , r ,<|> ) 


N 

(B cosh r - C sinh r) sin 0 + A n e nr sin n0 

n=2 


sin <t> 


J 


where N is a positive integer greater than two and A n , B, and C are arbi- 
trary constants. A two-dimensional transformation to the x,y plane is obtained 
by setting 0=0. 

Lines of constant r are transformed to circles in the x,y plane as r 
increases without limit in the negative direction. Because terms involving 
e nr vanish and -sinh r approaches cosh r as r increases without limit 
in the negative direction, equations (1) can be written approximately as 


x(0 , r) ~ C (B + C) cosh r] cos 0 


y(0,r) * C (B + C) cosh r] sin 0 


(r « 0) 

> 

(r « 0) 


( 2 ) 
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The line segment r = 0, 0 = 6^ 2V is transformed into a two-dimensional 
body in the x,y plane and becomes one boundary of the computational space in the 
0,r plane. An axisymmetric body is defined by mapping the line segment r = 0, 

0 = 0 = tt to the x,y plane and rotating the image around the x-axis 2tt radi- 
ans. Thus, 
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will therefore form two boundaries of the 
computational space for two-dimensional problems. The 0 = 0 plane, 0 = tt 
plane, 4> =□ 0 plane, and c{) = 2 tt plane are the boundaries in computational 
space for axisymmetric bodies. 


The final boundary of the computational space is chosen as r = r^. It 
will be shown in a later section that r^ can be mapped to negative infinity 
with one additional transformation. Such a transformation will cause the 
inflow-outflow boundary to be mapped to a circle of infinite radius as seen 
from equations (2). This transformation will simplify the problem of specifying 
inflow and outflow boundary conditions. An example of a typical transformed 
body shape is shown in figure 1 with the associated coordinate system. 


It should also be noted that the transformation can be written as a con- 
formal mapping for two dimensions. By setting 


w =□ x( 0 , r) + iy (0 , r) 
D = (B + C)/2 
A 1 = (B - C)/2 
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and simplifying, the following equation is obtained: 

N 

w = De iz - )> A n e - i nz 
n=1 

where z = 0 + ir and i ii\J-T 

The definition and orientation of the velocity components are determined 
by the coordinate system. The u component of velocity is directed along a 
line of constant r in the direction of increasing 0 . The v component of 
velocity is directed along a line of constant 0 in the direction of increas- 
ing r. Thus, from figure 1, it can be seen that the direction of positive u 

is counterclockwise around the body and the direction of positive v is toward 
the body. 

The detailed derivation of the metric coefficients for this coordinate 
system appears in appendix A. It is important to note that the metric coeffi- 
cient hg is equal to the metric coefficient h r . This equality greatly 

reduces the number of calculations that will have to be made when solving the 
governing flow equations. 


Determination of Body Shape 


Methods for determining a large variety of body shapes are discussed in 
detail in reference 19. For the sake of completeness, the methods for deter- 
mining the body shapes used herein are presented. 


For all A n = 0, the transformation from the 9,r plane to the x,y plane 
results in families of ellipses. For B = C, the transformation produces a 
circle of radius B. For B t C, the transformation results in an ellipse of 
eccentricity e c where 



Calculations have been made for configurations approximating planetary 
probe vehicles as shown in figures 2 and 3. The true body shape in figure 3 
is taken from reference 20. These approximations were obtained by considering 
the following parametric equations for a body shape 


4 

x b = C cos 0 - A n 
n=2 

4 

y b = B sin 0 + A n 
n=2 




cos n0 


> 

sin n0 


( 5 ) 


which correspond to equations (3) and (4) with N = 4 and <J> = 0. 
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In figure 4 some basic geometric parameters of a planetary probe are pre- 
sented. The radius of curvature of any point on the body can be expressed as 


R(0) 


r 2i 3/2 

|_1 + (dy/dx) J 


[Ox/99 ) 2 + Oy/99 ) 2 ] 3/2 

d 2 y/dx 2 


(9x/99)(92y/ae2) _ Oy/99)O2 x /302) 


( 6 ) 


By defining a nose radius of curvature %, a base radius of curvature R B , a 
vehicle length L, a maximum body radius Y mav , and the length from the nose to 
the location of the maximum body radius on the symmetry line Lq , it is possi- 
ble to solve for A 2 , A 3 , A 4 , B, and C. To specify Y mav , it is also nec- 
essary to determine 9 ma v - so that yt>(®max) = Y max and d Yb/ d ®(®max) = °* 
These specifications result in the following equations: 

(B + 4 A 4 + 3 A 3 + 2A 2 ) 2 

- ; = R(0) = R b (7) 

-C + I 6 A 4 + 9 A 3 + 4A 2 


( -B + 4 A 4 - 3 A 3 + 2A 2 ) 2 

— — - □ R(-tt) = R N (8 ) 

C + I 6 A 4 - 9 A 3 + 4A 2 

2C - 2A 3 = L (9) 

A 4 sin 49 max + A 3 sin 3®max + A 2 s ^ n 2 ®max + B s ^ n ®max = Y max (10) 

4 A 4 cos 49 raax + 3 A 3 cos 39max + 2A 2 003 2 ®max + B cos (®max) =9 (11) 

-A 4 cos 49 max - A 3 cos 3 9 max ~ A 2 coa 29 max + ^ ®max 

+ A 4 - A 3 + A 2 + C = L-| (12) 


Equations ( 7 ) to (12) can be solved by use of Newton's method to obtain 
A 4 , A 3 , A 2 , B, C, and 9 max . When choosing the values of R^, R B > a nd L-] , 
it may prove necessary to pick adjusted values which will give a better overall 
approximation to the desired body shape. As pointed out in reference 19, it is 
not possible to accurately model all body shapes with equations (3). Therefore, 
once the constants are obtained, the analytic approximation has to be compared 
with the desired body shape to determine whether agreement is satisfactory. A 
trial-and-error approach was established whereby values for Rjj, Rg, and Li 
were varied over a narrow range of values which were close to the measured val- 
ues for the desired body shape and then the resulting analytic approximation 
was compared with the true body shape. The analytic approximations shown in 
figures 2 and 3 were obtained in this manner. The parameters which were speci- 
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fied and the resulting values for B, C, A 2 , A 3 , and A 4 that correspond 

to these figures are presented in table I. 


ANALYTIC DEVELOPMENT 

In order to test the usefulness of the coordinate system just described, a 
solution to the supersonic flow at an angle of attack of 0 ° over blunt, axisym- 
metric, or two-dimensional bodies is obtained. The analysis is simplified by 
making the following assumptions: 


Np r = Constant 

p B = 0 

y = y(T) (using Sutherland's law) 

Perfect gas 

No radiative heat transfer 

The quantity yg is the bulk viscosity coefficient which is set equal to zero 
in keeping with general practices. 

Strictly speaking, the Navier-Stokes equations refer only to the momentum 
equations. By following Vincenti and Kruger (ref. 21), the term Navier-Stokes 
equations is used herein to refer to the complete set of partial differential 
equations that describe the motion of a viscous heat conducting fluid. These 
equations, expressed in a generalized orthogonal coordinate system, are nondi- 
mensionalized as follows: 

h = h*/R N « Y = Y*/R n * t = R N */V ro *t* 

p = p*/p * u = u*/V * v r v*/V * 

K K ^OO OO 00 

I = I*/V * 2 p = p*/p *V * 2 y = y*/y * 

The Navier-Stokes equations written in the transformed coordinate system 
may be obtained from references 22 and 23. Since only the steady-state solu- 
tion to the governing equations is desired, a term involving the time rate of 
change of pressure has been omitted from the energy equation in order to sim- 
plify the numerical procedure. The other time derivative terms must be kept in 
order to retain the hyperbolic nature of the governing equations in time. The 
omission of the transient pressure term prohibits an accurate modeling of the 
transient flow situation. However, as noted by Crocco (ref. 24), the unsteady 
equations are used strictly as a guide in choosing an iteration procedure. 

Allen (ref. 15) notes that the difference scheme can be modified to approximate 
a different unsteady equation as long as the steady-state difference equations 
and boundary conditions are correct. 
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The convective terms in the governing equations are written in conserva- 
tion form as recommended in reference 25 for shock capturing whereas the dissi- 
pative terms are expanded so that they no longer retain conservation form. 


In preliminary numerical calculations it was found that finite-difference 
approximations of terms involving products of metric coefficients with the 
conservation flow variables (that is, 9(hPu)/90, 9[ji 2 (p + Pu 2 )J/30, and 

9(h 2 Puv)/90) caused large errors in the numerical solution. This problem was 
most obvious in the free-stream flow ahead of the shock. Values for free- 
stream density along the stagnation streamline were in error by as much as 
4 percent. The corresponding Mach numbers differed from the free-stream con- 
dition by as much as 10 percent for = 2 and Y = 1.4. The governing equa- 
tions were rewritten in a form which separated derivatives of metric coeffi- 
cients from derivatives of conservation flow variables. The derivatives of the 
metric coefficients are analytic functions of 0 and r which are easily 
obtained from equations (A4). The derivatives of the conservation flow vari- 
ables are approximated by finite-difference formulas. The governing partial 
differential equations expanded in this manner become for: 

Continuity: 


3p 

3t 


-1 
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1 3n 
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(13) 


6 -momentum: 
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r-momentum: 
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Energy : 
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( 16 ) 


(Equation continued on next page) 
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The numerical solutions obtained by differencing the equations in this manner 
brought all free-stream properties within 1 percent of an undisturbed field. 

It should be noted that equations (13) to (16) are not in strict conserva- 
tion form in the sense that when they are applied to an undisturbed uniform 
flow, they do not return a uniform flow after one iteration. An error on the 
order of (A0)2 + (Ar)2 is introduced. For example, a uniform two-dimensional 
flow written in the coordinate system as expressed by equations (1) with B = 1, 
C = 1, A n = 0 has the following properties: 
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r 


1 



u = -sin 9 


v = -cos 6 
P = 1 
h = e -r 


Y = 1 


After one iteration (predictor step) on these conditions, the continuity equa- 
tion yields 


pn-1 = p" - Aterf =2H > (6 * A6 > ^ 3l 

pn+l n p n + Ate r cos q| ' sj ' n ^ — - 
pn+1 _ pn + At e r cos 0 


sin (0 - A0) 


cos 0 


1a0 2 A0 4 


3! 5! '7 

These errors have shown no sign of seriously disturbing the final solution. 


An additional coordinate transformation is utilized which simplifies the 
treatment of the inflow and outflow boundary conditions and gives some control 
on the density of mesh points near the body in a direction normal to the body. 

A new coordinate n is defined so that 

r = 3 log e n (n = e r/ £; 3 > 0) (17) 

where -°° < r ^ 0 and 0 < n = 1 . This coordinate stretching maps the inflow 
and outflow boundary conditions to infinity. All boundary conditions at infin- 
ity are known. This coordinate stretching does not affect the orthogonality of 
the coordinate system. Such a transformation was used by Kitchens (ref. 17) in 
order to eliminate rarefaction waves in a wake induced by fixing free-stream 
boundary conditions at a finite distance downstream of a body. Derivatives 
with respect to r are written as derivatives with respect to rt as follows: 
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All derivatives of flow properties with respect to r in equations (13) to (16) 
are now replaced by their counterparts in equations (18). 

On the line of symmetry the 0-momentum equation is not solved because 
u = 0. For axisymmetric problems, the limiting forms of the conservation equa- 
tions for mass, r-momentum, and energy must be obtained because the factor 1/Y 
which appears in many of the terms in those equations increases without limit 
as 0 approaches 0 or tt. (Note that Y = 1 for two-dimensional problems.) 

The limiting form of these equations is easily obtained by applying L’ Hospitals 
rule to all terms in these equations which are multiplied by 1/Y. 


NUMERICAL TECHNIQUE 

A modification of the Brailovskaya scheme which was introduced by Allen 
and Cheng (refs. 15 and 26) was applied to equations (13) to (16) to obtain the 
numerical solutions presented herein. The differencing technique is presented 
in appendix B. This particular method was chosen because the viscous stability 
limit in the Navier-Stokes equations is removed in Cartesian coordinates for 
constant viscosity. Allen and Cheng point out that when there is a large expan- 
sion of the flow around a corner or when the density in the near wake is small, 
the stability limit on the time step based on the diffusion terms can be 
severely restrictive as compared with the stability limit based on the inviscid 
terms. The elimination of this stability limit makes their modification to 
Brailovskaya’s scheme very attractive. The stability limit on the time step 
for the inviscid portion of the governing equations is derived in reference 27. 
This limit can be written in Cartesian coordinates as 

1 

At S - (19) 



In the transformed coordinate system, this limit is written as 

1 

At r — (20) 
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This transformed stability limit was obtained by substituting the terms 
for velocity and arc length in the transformed plane that correspond to the 
equivalent terms in the Cartesian plane. The stability limit from equa- 
tion (20) was used in all calculations without any indication of numerical 
instability. Also, comparisons were made in regions of severe expansion 
(around the corner of the Viking aeroshell, for example) and it was found that 
the actual time increment used exceeded the viscous limit by more than a factor 
of 20 (see table II), where the viscous stability limit is (ref. 27) 
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Convergence Criteria 

The modified form of the Brailovskaya method is used herein to obtain the 
steady-state solution of the governing equations. Because the solution is 
approached in an iterative manner, it is necessary to establish some criteria 
to insure that the results have in fact converged to a steady state. A vari- 
able e is defined so that 


e 



max 

i, j 


( 22 ) 


The variation of e as a function of iteration number is shown in fig- 
ure 5 for different mesh sizes. The solution is said to be converged when the 
value of e becomes less than some small number e. Care must be exercised in 
choosing the value of e as can be seen in figure 5. The value of e should 
be chosen as some small number which is below the level where the e variation 
begins to level out. For example, in the case of the e variation for a grid 
of 51 x 50 with the stretching parameter $ = 1, it can be seen that e starts 
to level out after 2000 iterations at 0.0002. When it was finally established 
that the variation of e had leveled and was starting a much slower decrease, 
a value of e equal to 0.0001 was chosen as the convergence criteria. A com- 

3i 

parison of results for — is given in table III at various iteration levels 

b 

to show the effect that various values of e have on the static enthalpy 
derivative at the wall. The enthalpy derivative should be very sensitive to 
any changes occurring in the flow field. The results from table III indicate 
that for a mesh size of 51 x 100, there is less than a 2-percent difference 
between the solution after 6000 iterations and the solution after 4000 
iterations. 


Fourth-Order Smoothing 

Preliminary calculations indicated that some kind of damping or smoothing 
routine was necessary to eliminate numerical instabilities, especially in the 
vicinity of the shock. A nonphysical damping function was used to eliminate 
these instabilities. (See ref. 1.) Terms of fourth order in the spatial grid 
are used to smooth results after every iteration as indicated in the following 
equations . 
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( 23 ) 




The tilde symbol over the flow variables indicates the undamped results from the 
second step of the difference scheme. 

The five-point formulas for the fourth-order derivatives have the form 


3 4 ( P n ) 


~as~ 


= p 1 ? 


(24) 


When the five-point formula for the fourth-order derivative involves property 
values across a line of symmetry, then the appropriate symmetric or antisym- 
metric value of the property must be substituted into equation (24); otherwise, 
nonzero values for 3p/3G and 9p/39 on the symmetry line will result. 


For j = 1 or j = NJ, no damping is used in the n-direction. For j = 2 

3 4 ( ) 

cannot be used because 

an 4 

it involves points outside of the computational space. Four-point formulas for 
third-order derivatives in the n-direction are used when j = 2 or j = NJ - 1 . 
These formulas have the form 


or j 2 NJ - 1 , the five-point formula for (An) 4 
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in equations (23) ; 


from equation (25) is substituted 
when j = NJ - 1 , the expression 
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in equations (23). 


from equation (26) is substituted for the term 
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In a stability analysis on the inviscid equations, Barnwell (ref. 1) shows 
that the damping coefficient must satisfy the inequality 


0 £ e S 1/24 


(27) 


Ideally, the smallest value of £ in this range which allows a stable smooth 
solution should be used. However, because of the large amount of computer time 
required to run a problem to convergence, no systematic search for an "ideal" 
value of £ was undertaken. It was found that £ = 0.001 would result in 
stable solutions for most of the cases considered herein. 


Initial Conditions 


Transformati on of unif orm velocity field.- In order to establish realistic 
initial conditions, it is necessary to resolve a uniform flow in the x,y plane 
to the equivalent condition in the 9,r plane. This can be accomplished in the 
following manner. Let a uniform velocity field of magnitude V m approach a 
body in the x,y plane at 0° angle of attack. All angles are measured from 
6=0 in a counterclockwise direction. Let be defined as the angle of the 
vector tangent to a line of constant r in the direction of increasing 0. 

From consideration of figure 6, it may be determined that 


u = V m cos 


v = -V^ sin 


( 28 ) 


It is now necessary to determine cos and sin \p as a function of 9 
and r. 


The differentials dx and dy along a line of constant r can be deter- 
mined from equations (1) as 


dx = 


N 


dy = 


(B sinh r - C cosh r) sin 0 + nA n e nr sin n0 

n=2 


N 


d0 


( 29 ) 


(B cosh r - C sinh r) cos 9 + nA n e nr cos n0 

n=2 
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From figure 6 and equations (29) and (A4), it can be shown that 


cos 


* = 


dx 


vdx 2 + dy 2 


N 

(B sinh r - C cosh r) sin 9 + nA n e nr sin n0 

n=2 
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N 


\ (30) 


sin = 


ay 


(B cosh r - C sinh r) cos 0 + y* nA n e nr cos n0 

n=2 


^/dx 2 + dy 2 


Thus, the initial uniform flow field may be resolved into its components u 
and v in the transformed system by using equations (28) and (30). 

Initia lization procedure .- When the first set of calculations were per- 
formed on the inviscid equations, the initial conditions were prescribed as a 
uniform velocity field everywhere except at the body where v was set equal 
to 0. This situation corresponded to a body suddenly materializing at rest in 
a supersonic stream at time t = 0. Problems developed as the shock started to 
form near the body. Negative values of static enthalpy occurred between the 
shock and the forebody; thus, any further calculations are terminated. A prob- 
able cause of these negative values was in the way surface boundary conditions 
are obtained. In most "shock capturing" numerical methods, in which no dis- 
crete shocks are assumed but rather are allowed to be smeared over several mesh 
points, it is common to see a slight overshoot and undershoot of properties on 
either side of the shock (ref. 25). An extrapolation procedure was used to 
obtain pressure on the body but this extrapolation was invalid when the proper- 
ties used in the extrapolation were in the region of the developing shock. 


Thus, either the method of obtaining the boundary condition had to be 
changed or the initial flow conditions had to be altered. Kitchens (ref. 17) 
took the first approach and allowed the surface velocity components to vary 
linearly from the uniform flow values to zero after some period of time t^. 
This specification corresponds to a porous body with the amount of mass flux 
through the surface decreasing to zero between times t = 0 and t = t-|. The 
no-slip condition is then specified on the surface for times t ^ t^. However, 
it was found to be simpler in the present program to alter the initial condi- 
tions in such a way that the shock would develop further away from the body. 
This was accomplished by calculating the transformed velocity components and 
multiplying the result by 1 - Then, the total enthalpy at a point was set 

equal to the free-stream static enthalpy plus one-half of the velocity squared. 
Thus, 


u = (1 - n )V co cos 
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The variable n = 1 corresponds to the body so that the no-slip condition 

is imposed. When q = 0, the velocities are unchanged from the uniform flow 

conditions. In the 6,p computational plane the initialized velocities vary 
linearly from p = 0 to p = 1. In the x,y plane this alteration corresponds 

to an exponential damping of the uniform velocity field as the body is 

approached from the far field. This specification does not model any physi- 
cally realistic situation but does allow a steady-state solution to be 
obtained. 


Boundary Conditions and Mesh Positioning 

An example of positioning of mesh points in the 0,p computational plane 
is shown in figure 7. Note that mesh points are positioned on the line seg- 
ments 0 = 0 and © T 3 and that mesh points are positioned one-half cell 

width away from the body, p = 1, and the inflow and outflow boundary, p = 0. 
This is the same technique used by Allen (ref. 15) and some explanation for the 
reasoning behind the approach will be given. Allen begins the derivation of 
the numerical approach by writing the conservation laws of mass, momentum, and 
energy in integral form. He works directly from the integral form of the equa- 
tions to obtain all the difference approximations. It should be noted that the 
same set of difference equations presented herein could be obtained from the 
conservation laws of mass, momentum, and energy written in differential form. 
However, the integral approach gives an additional insight to the derivation of 
the difference equations which can be lost if various second-order accurate 
difference approximations to derivatives that are obtained from Taylor series 
expansions are routinely substituted into the differential form of the conser- 
vation equations . 

All values of properties at a mesh point represent the average value of 
the property in the cell. Flux terms across the sides of a cell are written in 
terms of the average values. For example, 


Pu i 
i+- 
2 


» J 


= rl (P u )i+i , j + (p u ) 



where pu represents the average value in the cell and pu is the mass flux 
from the i cell to the i + 1 cell. 

When these values for fluxes and derivatives at the boundaries of a cell 
are substituted into the integral form of the conservation laws, the difference 
formulations of the conservation equations are obtained. For example, when sub- 
stituting into the integral form of the equation for the conservation of mass 
in Cartesian coordinates, one finally obtains 
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As can be seen, a second-order accurate central difference approximation 
results for the partial derivatives 9pU/3x and 3pV/3y which, as stated 
before, could have been obtained directly from the differential form of the 
conservation laws. In like manner the second derivative of U with respect 
to x would be written 


/a 2 u\ 

_ A 


. 1 . 
1+2 , J 

1 

Ui+1,j " U i,J 

U i,j - 

[d^J. . 

x / i, j 

1 X 
1 < 

it 

( Sx / 

. 1 . 
i~2» J 

= Ax 

Ax 

Ax 


- 2°i,J + 
Ax 2 


which is the second-order accurate central difference approximation to that 
derivative. 

In deriving the difference approximations to the derivatives in the gov- 
erning equations, the values of fluxes and derivatives -of velocities on the 
sides of the computational cell are the basis of the approximation. Thus, 
Allen (ref. 15) points out "The greatest advantage of the integral formulation 
is the conceptual aid it gives in applying boundary conditions. It indicates 
that wall boundaries would be better placed along cell edges rather than 
through mesh points . " 
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Consequently, in the generalized orthogonal coordinate system used herein 
all the property derivatives which appear in equations (13) to (16), except for 
6 derivatives on a line of symmetry, were approximated on the basis of the 
values of properties on the sides of a cell wall. All derivatives of metric 
coefficients are evaluated analytically. In most cases, that is, for all the 
interior points, the values at the cell wall are calculated by using the aver- 
age values of properties in the cells on either side of the wall. Methods for 
establishing cell wall values at the boundaries are discussed in appendix C. 


Symmetry Conditions 

The computational boundaries along the lines 0 = 0 and 9 = it corre- 
spond to a line of symmetry in the flow. The cell centers are placed on these 
boundaries to facilitate the calculations of the limiting form of the governing 
equations for axisymmetric flow situations. The symmetry conditions are 
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calculations it was 

found that 


line solution with smoother property distributions along the body could be 
obtained if the derivatives with respect to 0 were calculated by using a 
fourth-order accurate central difference formula. If a second-order accurate 
central difference formula is used to obtain derivatives with respect to 0 
across the symmetry line, then 9(pu)/90, 9(puv)/90, and 3u/39 are calcu- 

lated by using information from only one adjacent cell. A fourth-order accu- 
rate central difference approximation to the 0 derivatives feeds information 
from two adjacent cells to the symmetry line calculations. At i = 1, for 
example, 
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This differencing scheme does not follow the integral approach explained in the 
earlier sections. However, it was felt that the improved solutions obtained in 
this manner justified the special treatment of the symmetry line. 



RESULTS AND DISCUSSION 


A computational technique based on the numerical method described in the 
previous section has been programmed for use in the transformed coordinate sys- 
tem. The program has been written to describe supersonic flow over two- 
dimensional or axisymmetric body shapes. Comparisons with experimental data 
are presented which verify the technique. Flow fields around the Viking aero- 
shell and a candidate configuration for a Jupiter probe are also calculated and 
presented. A discussion of the program’s range of applicability along with the 
presentation of some computational results follows. 


Range of Applicability 


Practical limitations on storage space and execution time on the computer 
restrict the number of mesh points which can be used to define a flow field. 
The program requires 75 000 octal words of central memory for a mesh size of 
51 x 50 and 121 000 octal words of central memory for a mesh size of 51 x 100. 
Execution time on the Cyber 175 with an optimization for fast execution is 
approximately 0.9 sec per iteration per 1000 grid points. 


The stretching parameter 8 makes it possible to concentrate many grid 
points near the body at the expense of producing a sparse distribution of mesh 
points far away from the body. Some examples of the effects of the stretching 
parameter on mesh positioning are presented in table IV. For a fixed number of 
mesh points in the ri-direction , it can be seen that as 8 is decreased, the 
outermost mesh points are drawn closer to the body. Consequently, a lower 
limit on 8 must be imposed in order that the shock may be captured and the 
flow in the wake computed. 


By using a Blasius type solution for cases where boundary-layer theory 
applies, it can be shown that the nondimensional laminar boundary-layer thick- 
ness over an axisymmetric or two-dimensional blunt body is on the order of 



An adequate resolution of the boundary layer requires the 


placement of several mesh points across the actual boundary-layer thickness. 
Therefore, a large number of mesh points in the n-direction and a small value 
of the stretching parameter 8 are needed to position enough grid points in 
the boundary layer for very large Reynolds number flows. As indicated earlier 
in table IV, the use of coordinate stretching alone may pull all the mesh 
points too close to the body and prevent an adequate description of the bow 
shock shape and flow in the wake. For example, some typical flow field results 
indicate that the bow shock is stretched over three to four grid points. 
Morduchow and Libby (ref. 28) show that if one assumes that continuum assump- 
tions apply in a shock that is no less than seven mean free paths thick, then 
the Mach number of the flow must be less than 1.3- They express the mean free 
path of a molecule ahead of the shock as 
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By expressing 


Poo* N Re,oo 

the Knudsen number Nj^ n is obtained as 
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Since the nondimensional shock thickness is on the order of the Knudsen number, 
this relationship enables one to make an approximation of the size of the 
region over which the shock is smeared. Admittedly, the governing equations 
used herein are not sufficient to resolve the flow through the shock. However, 
resolution of the bow shock is needed to avoid smearing the shock throughout a 
large portion of the shock layer. Therefore, for the purposes of resolving the 
flow details with the tools of this analysis, computations are restricted to 
low Reynolds number flows. 


The free-stream Reynolds number, however, cannot be allowed to become too 
small or continuum theory, upon which the Navier-Stokes equations are based, is 
not valid. Probstein defines a Knudsen number in terms of a mean free path in 
the shock layer divided by the shock-layer thickness, 
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and uses this parameter to roughly delimit various flow regimes (refs. 29 
and 30). In particular, he defines two flow regimes which are within the capa- 
bilities of this analysis and place a lower limit on the Reynolds number range. 
In the viscous layer regime "the shock layer is a fully viscous continuum ame- 
nable to treatment with the complete Navier-Stokes equations , and the shock 
wave may be treated as a discontinuity across which the Hugoniot (shock) rela- 
tions apply" (ref. 30). This regime is delimited by 
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In the incipient merged layer regime, the shock layer is still treated as a con- 
tinuum but the shock can no longer be considered a discontinuity satisfying the 
classical Rankine-Hugoniot relations. This regime is delimited by 
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On the other hand, when 
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the fully merged layer regime is entered and it is incorrect to assume a con- 
tinuum state although Probstein cites several papers which indicate that the 
Navier-Stokes equations can be used to obtain estimates for skin friction and 
heat transfer, and to obtain mean density and velocity distributions. 

It should be noted that as the flow expands around the body into the wake, 
the mean free path can increase by approximately a factor of 20; this increase 
corresponds to the change of y*/p*\|T* from the forebody stagnation region to 
the wake. Thus, one could be in the viscous layer or incipient merged layer 
regime in the forebody and enter the equivalent of a fully merged layer regime 
in the wake. 

In order to keep the present analysis simple, the no-slip wall boundary 
conditions on velocity and temperature have been applied. In fact, the condi- 
tion for zero velocity slip has been used to establish the specialized form 
of the normal velocity derivatives in the vicinity of the wall. Probstein 
(ref. 30) derives a condition for negligible slip as 
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In all the cases considered herein, T s */Tb* ^ 1. However, it is empha- 
sized that this relationship does not mean that slip does not exist as an 
identifiable phenomenon but only that its effect on mean aerodynamic quantities 
is not large. Schaaf and Chambre (ref. 31) roughly bound a slip flow regime on 
the basis of experimental evidence by 
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Some of the experimental data used as a basis of comparison in the next section 
fall within this range and modifications to the solution technique are dis- 
cussed in that section. A more detailed examination of the effects of slip in 
this Mach number and Reynolds number range will be presented with the results 
of that comparison. 

A Reynolds number range may therefore be summarized as follows: 

0(100) S N Re>QO S 0(1000) 
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where the upper limit is imposed because of an inability to resolve a thin 
boundary layer and the lower limit is imposed so that continuum theory is 
applicable. A correction for velocity and temperature slip may be required 
at the lower Reynolds number limit, depending on other properties of the flow 
field. 

The process of capturing a strong shock and spreading the discontinuity 
over several mesh points can result in a small overshoot and undershoot of 
properties on either side of the shock. The magnitude of the overshoot 
decreases with increasing damping factor and increases with increasing mesh 
point separation. The bow shock is closest to the body, in a region of rela- 
tively dense mesh spacing, in the stagnation region; as the bow shock wraps 
around the body and approaches the Mach angle, it moves away from the body 
into an area of relatively sparse concentration of mesh points . It has been 
observed that as the bow shock forms in high Mach number flow, negative static 
enthalpies are calculated in the vicinity of the shock. These negative static 
enthalpies are a consequence of the velocity overshoot on the inflow side of 
the shock . 

The static enthalpy is computed by subtracting (u 2 + v 2 )/2 from the 
total enthalpy I. The total enthalpy is approximately constant across the 
shock and is equal to 


1 1 
= (Y - 1 )M 00 2 + 2 

By defining V-| as the calculated velocity ahead of the shock, the per- 
cent overshoot of the velocity Op is defined as 

(Vi* - V TO ») 

0 n = 100 = 100 (Vt - 1) 

p v * 1 

oo 


The calculated static enthalpy ahead of the shock is then expressed as 


il 



1 

(Y - 1 )M 00 2 
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2\100 + j (y . -| )M oq 2 100 20 000 


Consequently, if 


100 

Op > ~ 

(Y - DM* 2 

a negative static enthalpy is calculated and the program execution is 
terminated. 

As can be seen in these relationships, the allowable percent overshoot 
decreases as M^ increases. It has been found that with a grid size of 
51 x 50, the onset of the problem of calculating negative static enthalpies 
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ahead of a bow shock occurs at M^ ~ 4 for a sphere and 3*5 for a 

cylinder with y = 1.4. The Mach number limit is slightly higher for a sphere 
because the bow shock stands closer to the body in an area with a relatively 
large concentration of mesh points. 

It should be noted that improvements to the numerical approach can be made 
to extend the Mach number and Reynolds number range and to reduce the execution 
time on the computer. The technique can be reprogrammed to run on the STAR 100 
computer with a potential execution time reduction of an order of magnitude 
(ref. 32). The numerical approach can also be altered so that the stability 
limit on At is increased by using an alternating direction implicit method or 
partial implicitization (refs. 33 and 34). A shock fitting technique can be 
applied so that the bow shock is treated as a discontinuity which floats 
between grid points (ref. 35). This modification will extend the Mach number 
range and will remove the requirement of a large number of mesh points in the 
vicinity of the shock. Consequently, a more efficient coordinate stretching 
routine (ref. 36) can be used to include more points near the body and to com- 
pute higher Reynolds number flows. 


Experimental Comparisons 

Experimental pressure coefficient and convective heat-transfer coefficient 
distributions were obtained over a circular cylinder normal to a supersonic 
rarefied airstreara in the Mach number range 1.3 to 5 -9 > the Reynolds number 
range 19 to 2050, and at two cylinder wall average temperature levels of 90 K 
and 210 K in the low-density wind-tunnel facility of the University of 
California at Berkeley by Tewfik and Giedt (ref. 37)- Detailed comparison 
cases were run with the program described herein for a Mach number of 1.90 and 
a Reynolds number of 1 05 - 


Pressure coefficient comparisons for the case of an adiabatic wall are 
presented in figure 8. Close agreement is obtained between the calculated and 
experimental results over the forebody of the cylinder through a value of 

7T 

7T ^ 9 > — . The calculated pressure coefficients on the base for the range 
3 

TT 

0^9^- fall slightly below the experimental results but the comparisons 
3 

are still generally good here. The pressure coefficient over the cylinder for a 
cold wall case (T^ ~ 90 K) follows much the same trend as the adiabatic compar- 
ison (fig. 9). The densities near the wall in the wake for this case, and con- 
sequently the molecular mean free path, have changed by approximately a factor 
of 10 from the stagnation region to the base. 


A recirculation region is formed in the rear stagnation area of the cylin- 
der. The u component of velocity at the mesh point closest to the body 


changed sign at I 



The recirculation region is roughly bounded 


27 


by two stagnation points along the axis of symmetry in the wake, one at the 
body (x = 1 , y = 0) and the other approximately one body diameter downstream 
at x = 2.2, y = 0. A study of the onset of separation from the body in the 
base region is presented in appendix D. 

Convective heat-transfer coefficients are compared in figure 10 for an 
average wall temperature of 90 K and for a linearly varying wall temperature 
which more closely approximates the actual experimental temperature variation. 
These heat-transfer coefficients were evaluated as follows: 


Si/3n | b ~ C P* M * 1 3i 

Np r %*hB i aw - ib 


approximately 25 percent greater 
Calculations were performed with 
twice the number of mesh points in the normal direction, and with a stretching 
parameter of 0.5 to bring the mesh points even closer to the body. Changes in 
heat transfer were insignificant with the finer mesh (less than 4 percent of 
the coarser mesh solution). The results of figure 10 were used to compute con- 
vective heat transfer so that the actual heat-transfer errors could be measured 
The convective heat-transfer results are presented in figures 11 and 12 where 

3c, exp ° h c,exp( T aw,exp ~ T* , exp) 

Qc,calc = ^c , calc (^aw, calc ~ T b,calc) 

However, comparisons in the stagnation region are still in error by approx 
imately 25 percent and remain at about that level over most of the body. For 
low temperatures (=100 K) and at 1 atmosphere pressure, the National Bureau of 
Standards tables (ref. 38) show the Prandtl number at approximately 0.77 for 
air. The convective heat transfer for a linear temperature variation with 
Np r = 0.77 was then computed and the results appear in figure 12. Errors in 
the stagnation region have been reduced to approximately 18 percent but the 
overall comparison is still poor. The viscosities computed by Sutherland's 
equation were found to be within 5 percent of the experimental data (ref. 39) 
at the wall temperature range 80 K < T* < 110 K and the effects of e on the 
enthalpy derivatives were found to be small (see appendix E) ; therefore, these 
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As indicated in the figure, the results are 
than experimental results in the forebody. 


28 



factors could not account for the errors of 18 percent in the stagnation 
region. Consequently, it was decided that slip conditions must be important 
in this Mach number and Reynolds number range. 

A condition for temperature slip as a boundary condition was incorporated 
into the program as follows. The condition for temperature slip (ref. 31), 
velocity slip being ignored, can be expressed as 


T 0 * “ T b » = 


2 - a t 2y £* 9T* 

a b y + 1 Np r 9n* 


Jl« = 


U26f 



* 


(34) 

(35) 


In equation (34), is the thermal accommodation coefficient and T 0 

is the slip temperature. A value for a b = 0.9 has been found to be repre- 
sentative of values for thermal accommodation coefficients as presented in ref- 
erence 31. Equation (34) can be simplified by using equation (35), Sutherland’s 
law, and the equation of state to obtain 
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For zero velocity slip the static enthalpies in equation (36) can be replaced 
by the total enthalpies. By substituting the finite- difference formulation for 

ail 

from equation (C4), the slip value of wall enthalpy may be obtained 
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Convective heat transfer was calculated by use of this temperature jump 
condition and good comparisons were observed with the experimental data in fig- 
ures 11 and 12. For this Mach number and Reynolds number combination, it is 
clear that the proper slip boundary conditions are necessary for calculating 
heat transfer. 

Berezkin, et al. (ref. 40) applied holographic interferometry to the 
supersonic low density gas flow about a spherical body at a Reynolds number 
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of approximately 1750. Shadow and interference images of the density field 
were obtained in the shock layer of a spherical projectile at Mach 1.93 using 
a Mach-Zehnder interferometer. Density distributions are presented along two 
vertical lines corresponding to x = -1 , 0 < y < y s and x = 0, 1 < y < y s * 

Experimental values for shock standoff distance were not presented and all the 
distributions are nondimensionalized by the distance, y s - y^. The analytic 
value of shock standoff distance has been computed from the midpoint of the 
captured shock. Experimental and analytical results for the density variation 
are presented in figures 13 and 14. Agreement is good across most of the shock 
layer. Discrepancies near the wall occur because the wall temperature in the 
experimental case is not at a steady state. (An adiabatic wall was assumed for 
the purpose of making these calculations.) By imposing a constant wall temper- 
ature corresponding to = 1^ = 0.61^, the density agreed well near the wall 
at x = 0 but was too high near the stagnation point at x = -1 . No further 
attempts were made to approximate the unsteady nature of the wall temperature. 

The characteristics of the solution in the vicinity of a shock are well 
demonstrated in figures 13 and 14. The locations of the mesh points used to 
calculate these distributions are marked with crosses. These solutions were 
obtained with two different damping constants, e = 0.005 and e = 0.04. The 
results for the smaller value of e (fig. 14) show a sharper shock profile 
with more of a density undershoot as compared with the e = 0.04 results. 
Because of the nature of the coordinate system, the density profile along the 
line x = -1 had to be obtained by linearly interpolating between mesh points 
in the n-direction. The cross marks which appear on the shock in figure 13 
refer to interpolated mesh points and consequently distort the effect of e in 
that figure. 

An attempt has been made to verify a shock standoff distance calculation 
by comparing it with the experimental data for supersonic flow over a sphere 
with M ro = 4.2, NR e>00 = 200, T^* = = 300 K in nitrogen. An electron 

beam X-ray technique was used to obtain measurements of the density along the 
stagnation streamline of a spherical model (ref. 41). This Mach number is on 
the border of conditions for which the occurrence of a velocity overshoot can 
cause the calculation of a negative static enthalpy. In this particular case, 
for a mesh size of 51 x 50 and B = 0.5, it was determined on the basis of 
experimental evidence that there were approximately four mesh points across the 
shock on the stagnation line. Consequently, there were no problems in evaluat- 
ing static enthalpies ahead of the shock in the stagnation region. However, as 
the shock started to form around the body, it entered a region where mesh 
points became separated by larger distances. After approximately 350 itera- 
tions, a negative static enthalpy was calculated for this case at a value of 
2tt 

0 ^ — . A modification was made in the program logic which set the static 
3 

enthalpy equal to the free-stream static enthalpy whenever a negative value was 
computed. The density profile across the shock, after 2000 iterations with 
this modification, is presented in figure 15. Comparisons here show good 
agreement with experimental data. 
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Flow Fields Over Planetary Probes 


The analytic approximations which correspond to the Viking aeroshell and to 
a Jupiter probe as shown in figures 2 and 3 and defined in table I were used to 
compute the flow fields over typical entry bodies at = 2 and Np e m = 100. 
There are no experimental data currently available for these vehicles at this 
Mach number and Reynolds number range. The previous comparisons of calculated 
and experimental density distributions, heat-transfer distributions, and shock 
predictions indicate that the calculated flow field for the planetary probes 
should be representative of the actual fields. Flow fields were computed over 
the Viking aeroshell with a cold wall boundary condition, I5 = 0.61^, and an 
adiabatic wall boundary condition using y = 1.285, Np r □ 0.685, and 
T^* = 222 K. A flow field over the Jupiter probe was computed with an adia- 
batic wall condition and y = 1.667, Np r = 0.675, and T^* = 27*7 K. The grid 
size for all the cases was 51 x 50 with 8 = 1 and e = 0.04. 

A contour plot of density over the Viking aeroshell is shown in figure 16 
for the adiabatic case and in figure 17 for the cold wall case. Contour lines 
were obtained by using the program CONTOUR (ref. 42). CONTOUR prepares a 
Cartesian grid system based on interpolated values from the transformed coordi- 
nate system. Contour lines are then plotted based on interpolation from the 
Cartesian grid. Consequently, there is some distortion in the density contours 
due to the double interpolation and due to differences in the concentrations of 
mesh points between the Cartesian grid (which gives a uniform concentration of 
mesh points through the field) and the transformed coordinate system (which 
concentrates mesh points close to the body) . Irregularities in contour lines 
are caused by these distortions. 

The bow shock is shown clearly in both figures with the Mach angle 
approaching 60° as would be expected for M^ □ 2. Based on an argument that 
the nondimensional shock thickness is on the order of the Knudsen number, it 
appears that the calculated shock thickness is approximately 10 times the 
actual value. The shock in the stagnation region is smeared over three mesh 
points which is consistent with results obtained in the experimental comparison 
section. Because of the computational costs, no attempt was made to improve 
shock resolution by including more mesh points. The thermal boundary layer for 
the cold wall case (no temperature slip) is evident in the forebody in fig- 
ure 17. There are approximately five mesh points through the thermal boundary 
layer in this case. The Mach number and Reynolds number range for this case 
closely corresponds to the first comparison case in which good agreement with 
heat-transfer data was obtained by using the temperature slip condition. So, 
except for the effect of temperature slip, the thermal boundary layer should be 
accurately predicted. 

Velocity vectors over the Viking aeroshell are presented in figure 18 for 
the adiabatic case. The recirculation region is clearly visible in the wake. 
Here again, the effect of velocity slip on the recirculation region is subject 
to future verification. Pressure coefficients and heat-transfer distributions 
are presented in figures 19 and 20 for the cold wall case. The method used to 
integrate pressure coefficients and skin friction to obtain drag coefficients 
is presented in appendix F. The breakdown of the drag coefficient into the 
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forebody and afterbody pressure drag coefficient and the forebody and afterbody 
skin- friction drag coefficient is presented in table VI for the adiabatic case. 

Similar presentations are made for flow fields obtained, an adiabatic wall 
being assumed, around a Jupiter probe. Density contours are presented in fig- 
ure 21. Here again, the calculated shock thickness is probably on the order of 
10 times the true value. Velocity vectors are presented in figure 22. A mag- 
nification of this vector field with more velocity vectors included shows a 
distinct velocity gradient near the body from which a boundary-layer growth can 
be observed (fig. 22). Such a clear, boundary-layer type phenomena will cer- 
tainly become less obvious with the inclusion of velocity slip. A pressure 
coefficient distribution is presented in figure 23 and a drag coefficient 
breakdown is presented in table VI. 

When calculating the flow field over the Jupiter probe with stretching 
parameter 3 = 0.5, it was found that there were not enough points to define 
the wake flow far from the body. The last computational points in the wake 
before the outflow boundary were still in a region of subsonic flow. The tech- 
nique of mapping the outflow boundary to infinity when these last computational 
points were still in a subsonic region initiated large oscillations in the flow 
variables . The problem was eliminated with 3=1. 


CONCLUDING REMARKS 

The coordinate transformation procedure and solution technique described 
herein has been shown to be well suited for the calculation of the complete 
flow field surrounding various two-dimensional and axisymmetric bodies. This 
conclusion is supported by the good comparisons obtained between predicted val- 
ues and experimental data for pressure coefficient distribution and heat- 
transfer distribution on a cylinder and for density distributions and shock 
standoff distance on spheres. The solutions for the flow over the Viking aero- 
shell and a Jupiter probe were obtained by simply adjusting five transformation 
constants without having to perform any adjustment on boundary conditions. The 
versatility of the technique in being able to model many different axisymmetric 
blunt body shapes and its ability to calculate both the forebody and wake flow 
over these shapes represent a significant advance in aerothermodynamic 
technology. 

Improvements to the numerical approach can be made to greatly extend the 
Mach number and Reynolds number range and to reduce the execution time on the 
computer. The technique can be reprogrammed to run on the STAR 100 computer 
with a potential execution time reduction of a factor on the order of 30. The 
numerical approach can also be altered so that the stability limit on At is 
increased by using an alternating direction implicit method or partial 
implicitization . 

An upper limit on Reynolds number on the order of 1000 is imposed because 
of an inability to resolve a thin boundary layer with a reasonable number of 
computational mesh points. A lower limit on Reynolds number on the order of 
100 is imposed so that continuum theory is applicable. An upper limit on Mach 
number approximately equal to 4.0 is imposed because of a tendency to calculate 


32 



negative static enthalpies in the process of capturing the bow shock with a 
reasonable number of mesh points for large Mach numbers. These limitations 
are a consequence of the numerical method used and not of the coordinate system 
itself. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
October 26, 1977 
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APPENDIX A 


METRIC COEFFICIENTS 

The metric coefficients of an orthogonal curvilinear coordinate system are 
of the nature of scale factors which give the ratios of differential distances 
to the differentials of the coordinate parameters (ref. 43). The equations 
which are used to define the metric coefficients can be derived in the follow- 
ing manner. 

Consider the orthogonal transformation as defined in equations (1) for 
axisymmetric bodies. Let a 0 curve be defined in space as a line of con- 
stant r and <J> in equations (1). If ? is the position vector of a 
point in space, then a tangent vector to the 0 curve is given by 


^ 3r 3r 9s 0 

90 9sq 90 


(A1) 


where sq is arc length along the 0 curve. The vector 9r/9sQ is a unit 
vector tangent to the 0 curve which is redefined as uq (fig. 24). Equa- 
tion (AD can now be rewritten as 


$ = hQue 


(A2 ) 


ds 0 * 

where ha = is the length of 0. Therefore, the relation for he can be 

d0 

written 
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3x 3y -v 3z -> 

— i + — j + — k 

30 30 30 


(A3) 


The final expression for h0, along with h r and h^, which can be obtained in 
a similar manner, becomes 
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(Equations continued on next page) 
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Note that hg^ =_^h r ^ in equations (A4). The metric coefficients are all 
positive as long as^ U 9 is in the same direction as 0 and u r is in the 
same direction as R in a right-hand coordinate system. Therefore, the metric 
coefficients may be written as 


hg = h r = h 


(A5) 


h«j, = |y(0,r,O)| = Y 


(A6) 
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APPENDIX B 


MODIFIED BRAILOVSKAYA SCHEME IN TRANSFORMED COORDINATE SYSTEM 

The differencing procedure presented in reference (15) was applied to the 
governing equations in the transformed coordinate system. The subscripts 0 
and r in this section refer to derivatives with respect to 0 and r, 
respectively. Derivatives with respect to r are approximated by difference 
formulas in 0 as defined in equations (18). The difference approximations to 
the various derivatives and details of the iteration procedure follow: 
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The term C 2 is a constant which depends on the gas and the free-stream 
temperature. These constants can be found in reference 44. For example, 

C 2 = 100/TJ* (For air) 

C 2 = 233/T ra * (For carbon dioxide) 

C 2 = 98/T^* (For helium) 


The metric coefficients and their derivatives are all analytic functions of 0 
and r which can be obtained from equations (A4). The right-hand side of each 
conservation equation is now written as 
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The modified Brailovskaya scheme can now be written as 
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BOUNDARY CONDITIONS AND DIFFERENCE FORMS AT THE BOUNDARIES 


Inflow and Outflow Boundary 


As discussed in a previous section, the computational boundary at r> = 0 
corresponds to a circle of infinite radius surrounding the body. This one 
boundary has mass entering from the far field and leaving to the far field; 
hence, it is both an inflow and outflow boundary. Values for properties on the 
cell walls at this boundary are 

au ro 
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Values of derivatives with respect to n at mesh points whose cells lie 
along this boundary (j = 1) are calculated as follows: 
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In like manner, 
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Examples of second derivative formulas are 
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Wall Boundary Conditions 

The computational boundary at T) = 1 describes the entire body in the 
(x,y) plane. As a consequence of the no-slip condition u b = 0 on the cell 
walls that are bounded by h = 1 . Since there can be no mass flux through the 
body, v b = 0 on the cell walls that are bounded by ri = 1 . Values for pres- 
sure on the wall are extrapolated from interior points. Thus, 


Pb = 1*5Pi } NJ ~ 0-5Pi,NJ-1 


When the enthalpy or temperature of the wall is specified, the density at the 
wall can be calculated directly from the equation of state. For the case of an 
adiabatic wall, the enthalpy at the side of the cell is calculated from interior 


3l 

points at the previous time step using a backward difference formula for — 

3n 

which must equal zero for the adiabatic case. Thus, for an adiabatic wall, 
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nT n-1 T n-1 
9INJ " Inj-i 


L b" = 


(Cl) 


Mass, momentum, and enthalpy flux across the cell walls opposite to the 

1 

wall boundary, corresponding to NJ , are calculated differently than those 

2 

for other interior cell interfaces. It was observed in reference 15 that very 
low, and, in some cases, negative densities would occur in cells near the wall 
above the separation point on the base for small mesh sizes and low Reynolds 
numbers. Allen and Cheng traced this behavior to the way mass, momentum, and 

1 

energy flux were calculated along the NJ - — line of cell edges. They noted 


that v b a 0 and that as a consequence of the continuity equation and the no- 
3v [ 

0 in the steady state. This condition implies that 


slip condition, 


3n 


the variation of v as a function of n near the wall should be quadratic. 

1 

However, by defining the velocity and fluxes through the cell wall at NJ - — 

2 


as the average of values at NJ and NJ - 1 , a linear variation of velocity 
in that area is implied. It would be more accurate to determine the value of 

1 

velocities and fluxes at NJ - — on the basis of a quadratic variation near 


the wall by putting a second-degree polynomial through the three points at the 
two cell centers and at the wall. Consequently, 
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Values for the derivatives of u, v, and I with respect to r\ at the 
wall are obtained from a backward difference formula using interior points. 
Therefore , 
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Some examples of the difference approximations to derivatives with respect 
to n in cells positioned next to the wall (at j = NJ) are presented: 
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Because of the redefinition of flux terms through the cell walls at 
1 

NJ - — , the derivatives with respect to n at cell centers along the line 
NJ - 1 are defined as follows: 
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Other n derivatives are defined similarly at J = NJ - 1 . 
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OBSERVATIONS OF SEPARATION IN WAKE 

The presence of an inflection point in the pressure distribution curves 
can usually be used to determine the location of the boundary-layer separation 

3 2 p /Pi+1 - 2 Pi + Pi-1 

from a cylinder surface (ref. 45). Values for = I— 

30 2 \ A0 2 

were calculated in the vicinity of the separation point of spheres, cylinders, 
and the planetary probe shapes presented herein with cold wall and adiabatic 
wall specifications. In all cases it was found that separation occurred within 
+A0 of the location where 3 2 p/30 2 changed sign. It should be noted that 
this derivative is not equivalent to 3 2 p/3s 2 where s is the arc length. 
Since ds □ h d0 , it can be shown that 

3 2 p 1 /3 2 p 1 3h 

3s 2 h 2 y30 2 h 86 

3h 

where — =0 for a sphere or cylinder but is nonzero for any other body 
30 

shape. 
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EFFECT OF e ON ENTHALPY DERIVATIVE 

Because the convective heat transfer is so sensitive to the enthalpy 
derivative at the wall, it was decided to examine the effect of e on this 
derivative. Three flow field solutions were obtained for the first experi- 
mental comparison case (M^ tj 1.90, Np e = 105, Np r =0.7, y = 1.4, 

Tfc# = 90 K) . The first solution was obtained with e = 0.04, 8=1, and a 
mesh size of 51 x 50; the second solution, with e = 0.001, 8=1, and a mesh 
size of 51 x 50; and the third solution, with e = 0.001, 8 = 0.5, and a mesh 
size of 51 x 100. Results of this investigation are presented in table V. It 
can be seen that there is only an average of 2-percent difference between the 
e = 0.04 solution and the e = 0.001 solution for the mesh size of 51 x 50. 
Most of this error is due to differences in the derivatives on the rear stagna- 
tion line. However, the magnitude of these differences is small in the base 
and will not appreciably change heat-transfer levels there. The finer mesh 
size causes an average of 4-percent increase in the enthalpy derivatives . 

These derivatives are thus seen to be more sensitive to mesh size near the wall 
than to the smoothing parameter. 
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DRAG COEFFICIENTS 

The drag coefficient around an axisymmetric body in the transformed coor- 
dinate system can be expressed as 
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where the angle ^ is defined in figure 6 and sin ip and cos ^ are defined 
in equations (30). Also, 
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Combining equations (FI) to (F4) yields 
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By substituting the values of h and Y from equations (A4) to (A6) and 
the value of sin and cos ^ from equations (30) and noting that r = 0 for 

the body, one obtains 
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TABLE I.- PARAMETERS DEFINING PROBE SHAPES 


Parameter 

Viking 

aeroshell 

Jupiter 

probe 

R n 

1 

1 

r B 

2.5 

3.227272727 

L 

1.5454545 

3.227272727 

L 1 

Ymax 

.588940466 

1 .378281288 

1 .590909091 

2.163636364 

B 

1 .514647678 

2.078979297 

C 

.716078733 

1 .548980197 

a 2 

-.059786238 

-.042862760 

a 3 

-.056648540 

-.064656167 

A 4 

.042817883 

.062550354 


TABLE II.- VISCOUS STABILITY LIMIT CALCULATIONS 
ON VIKING AEROSHELL 
j^Atinviscid = 0.01100; ri - 1~j 




Atinviscid 

Index i 

Atviscid 

Atyiscid 

1 

5.575 X 10-3 

1.973 

2 

5.530 

1.989 

23 

.604 

18.212 

24 

.494 

22.267 

25 

.462 

23-809 

26 

.506 

21.739 

29 

1 .205 

9.129 

51 

35.703 

.308 


cr\ 
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TABLE III.- VARIATION OF ENTHALPY DERIVATIVE WITH INCREASING ITERATIONS 

[(51 x 100), e = O.OOl] 



1 

3i/3n, 

3i/9ri, 

3i/8ri, 

8i/3n, 

3i/3n, 

9i/3n, 

Index i 

0, deg 

6000 iterations, 

5000 iterations, 

4000 iterations, 

3000 iterations, 

2000 iterations, 

1000 iterations, 



e □ 8.45 x 10-5 

CD 

ii 

X 

o 

Jr 

e = 3.24 x 10" 4 

e = 6.19 x 10“ 4 

e = 1.40 x 10-3 

e = 4.42 x 10-3 

1 

0 

-0.653 

-0.660 

-0.660 

-0.707 

-0.900 

-1 .540 

6 

18 

-.580 

-.587 

-.587 

-.580 

-.707 

-1.307 

11 

36 

-.773 

-.773 

-.780 

-.787 

-.747 

-.813 

16 

54 

-1 .160 

-1.160 

-1.167 

-1.173 

-1.133 

-.760 

21 

72 

-2.027 

-2.033 

-2.033 

-1.987 

-1.953 

-1.467 

26 

90 

-3.533 

-3.533 

-3.480 

-3.493 

-3.407 

-2.780 

31 

108 

-5.640 

-5.640 

-5.647 

-5.600 

-5.513 

-4.853 

36 

126 

-8.293 

-8.293 

-8.300 

-8.307 

-8.207 

-7.547 

41 

144 

-10.993 

-10.993 

-10.993 

-10.987 

-11.000 

-10.367 

46 

162 

-12.987 

-12.980 

-13.033 

-13.080 

-13.073 

-12.527 

51 

180 

-13.687 

-13.747 

-13.733 

i 

-13.780 

i 

-13.827 

! 

-13.273 


i 

i 


mi mi i ii mi uni in 



TABLE IV.- EFFECTS OF B ON DISTANCE d OF MESH POINT 
FROM ORIGIN OF SPHERE OF RADIUS EQUAL TO ONE 


j 

3 = 1.00 

B = 0.75 

B = 0.50 

B = 0.25 


d 

d 

d 

d 

1 

100,000000000 

31,622776602 

10,000000000 

3, 162277660 

2 

53,533333333 

13,872638168 

5,773502692 

2,402811414 

3 - 

20,000000000 

9,457416090 

4,472135955“ 

2,114742527 

“T 

14,28571928b 

7.348123891 

3,779644730 

1,944130842 “ 

5 

11.111111111 

6,085806195 

3,333333333 

1.825741858 

6 

9,090909091 

5,235467719 

3,015113408- 

1.736408203 ” 

7 

' 7,692307692 

4.61894175S 

2,773500981 " 

1 ,665363.133 — 

8 

b, 666666667 

4, 1 4888651Y 

2,581988897 ” 

1 , 606856838 

”9 

5.882352941 

J, 777141971 

2,425356250 " 

1.557355531 

10 

5,263157895 

3,474839897 

2,294157339 

1,51 4547596 

U 

9, 761909762 

3.223558299 

2. 182178902 

1 .477219991“ 

12 

9,397826087 

3,010954016 

2,085144141 

1 ,4440028 19 

n 

4,OOOOUOOOO 

2.828427125 

2,000000000 

1,414213562 

”T4 

3,703703709 

2,669790460 

1,924500897 ' 

i ,387263817 

T5 " 

~ 3,448275802 

2.530471799 

1 ,B5695T362 “ 

1 73627irr768 — 

16 

3.2258-06452 

2,407014628 

1.796053020 

I ,340169027 

17 

3. 030303030 

2,296754337 

1 ,740/76560 

1 ,319384917 

18 

2,857192857 

2,197601621 

1,690308509 

1,300118652 

19 

2,702702703 

2,107893705 

1,643989873 

1, 262181665 

20 

2.564102S64 

2, 026289737 

1,601281538 

1,265417555 

21 

” 2,459024390 

“ 1,951695710 

1 ,56 1 7376 T9 ' 

1 ,2496950X0“ 

22 

2,325581395 

1 .883209595 

i, 524985703 

1,234905115 

23 

2,222222222 

1.820080575 

1 .490711985 

1,22094/167 

24 

£.127659574 

1.761678307 

i ,458649915 

1 ,207745799 

25 

2.040816327 

1,707469442 

1 ,4285 71429 

1,195228609 

26 

1 ,960784314 

1 ,656999464 

1 ,400260064 

1,183534306 

”27 

1,686792453 

“ 1,60-9878490 

I, 373605659“- 

I ,172009232” 

28 

1,818181818 

1.565770054 

1,348399725 

1,161 206 15 i 

29 

1,754365965 

1,524382162 

1.324532357 

1 , 150685294 

30 

1 .694915254 

1,485460101 

1 ,3018891 10 

1 ,141003554 

31 

l ,639344262 

1,448780606 

1 ,280368799 

1,131533826 

32 

1,587301567 

1,414147102 

1,259881577 

1 ,122446465 

“33 

1 ,538461538 

I .381385TB5 

I ,240547346 

1 , 113708825 

39 

1,492537313 

1,350342380 

1.221694444 

1,105302874 

35 

1 ,449275362 

1.320879445 

1.203858531 

1,0972040/2 

36 

1 ,408450704 

1,292874109 

1,186781658 

1,089395068 

37 

1,369863014 

1,266216169 

1.170411472 

1.081855569 

38 

1,333333333 

1.240806479 

1,154700538 ” 

1,074569932 

39 

T”, 2987 0 1 299 1,216555583 ' 

1 , 139605765 

1,067573192“ “ 

90 

1,265822785 

1,193382546 

1 • 125087901 

1,060701608 

91 “ 

I ,23456790 i 

1.171213948 

l.iiiiiiin 

1.054092553 

92 

1 .20*8*9277 

1.149983028 

1.097642600 

1,047684399 

93 

l t 176^70588 

1 . 129628929 

1,084652289 

1,041466413 

99 

1 ,139325267 

1 , 1 10096058 

1,072112535 

1,0354286/2 

95 — 

1,123595506 

l , 09 1 333521 - 

1.059997800 

1,029561985 

96 

1,098^01099 

1,073294628 

1,048284837 

1,023857821 

97 

1 , 0752688 1 7 

1.055936467 

1.036951695 

1,013308251 

98 

1,052631579 

1.039219521 

1,025978352 

1,012905895 

99 

1 ,030*27835 

1,023107337 

1,015346165 

1,007643868 

50 

1,010101010 

1.007566232 

1,005037815 

1,002515743 
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TABLE V.- EFFECT OF 


3 AND MESH SIZE ON ENTHALPY DERIVATIVE 



© 


© 


© 


e = 0.04 

e = 0.001 

e = 0. 


(51 x 50) 

(51 x 50) 

(51 x 

Index i 

3i 


3i 


8i 


in 

b 

in 

b 

in 

1 

-0.598 

- 0.688 

- 0.6 

6 

-.544 

-.566 

-.5 

11 

-.774 

-.783 

-.7 

16 

-1.139 

-1.134 

-1.1 

21 

-1.954 

-1.936 

-2.0 

26 

-3-324 

-3.283 

-3.5 

31 

-5.301 

-5.259 

-5.6 

36 

-7.809 

-7.799 

-8.2 

41 

-10.434 

-10.444 

-10.9 

46 

-12.453 

-12.499 

-12.9 

51 

-13.213 

-13.216 

-13.6 


100 ) 


Percent differences 


Average, percent 4.60 


1 nn 

® - © 

1 nn 

© 

1 1 

© 

1 uu 

© 

1 uu 

(D 

© 

-8.42 

5.36 

-13.08 

-6.21 

-2.41 

-3.89 

.13 

1 .29 

-1.15 

- 1 .81 

-2.24 

.44 

-3.60 

-4.49 

.93 

-5.92 

-7.08 

1 .25 

-6.01 

-6.76 

.80 

-5.84 

-5.96 

.13 

-5.09 

-4.99 

-.10 

-4.11 

-3.76 

-.37 

-3.46 

-3.44 

-.02 

, 4.60 

4.34 

2.015 


100 


TABLE VI.- DRAG COEFFICIENTS 



Pressure 

drag 

coe fficient 

Skin-friction 

drag 

coefficient 

Total pressure 
drag 

coefficient 

Total skin- 
friction drag 
coefficient 

C D 

(No slip) 


Forward 

Aft 

Forward 

Aft 


Jupiter probe 

1 .051 

0.1975 

0.1685 

0.06267 

1.243 

0.2285 

1.477 

Viking adiabatic wall 

1.296 

0.2610 

0.1121 

0.03650 

1.557 

0.1486 

1 .705 
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Figure 1.- Generalized natural coordinate system. 
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Figure 2.- Analytic approximation to Viking aeroshell with 
associated coordinate lines. 
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Analytic approximation 

Planetary probe vehicle 
(Candidate configuration for 
Jupiter probe, ref. 20 ) 




Figure k . - Basic geometric parameters of planetary probe. 




Figure 5.- Convergence test results for various mesh sizes 






























Figure 9-” Pressure coefficient distribution on cylinder with cold wall. 
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Figure 11.- Convective heat-transfer distribution on cylinder with constant wall temperature. 
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Figure 15.- Density distribution across shock layer 
on stagnation line with check for shock standoff 
distance. 
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s Arc length increases 
toward body 


Figure 2b. - Schematic of vectors which determine metric coefficients. 
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